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Abstract 


This  paper  describes  an  algorithm  and  a  FORTRAN  subprogram,  CHAIN, 
for  simulating  the  behavior  of  an  (n+1)  state  Markov  chain  using  a 
variance  reducing  technique  called  rotation  sampling.  The  simulation  of 
k  miaroreplioations  is  carried  out  in  parallel  at  a  mean  cost  ^0(ln  k) 

<v 

and  with  variances  of  sample  quantities  of  interest  £)0((ln  kj^/k^)  . 

The  program  allows  for  independent  maororeplioations  ,  each  of  k  micro¬ 
replications,  in  order  to  facilitate  estimation  of  the  variances  of  sample 
quantities  of  interest.  The  paper  describes  theoretical  results  that 
underlie  the  algorithm  and  program  in  Section  1  and  presents  applications 
of  interest  for  first  passage  time  and  steady-state  distributions  in 
Section  2.  Section  3  describes  the  algorithm  and  CHAIN  and  an  example 
in  Section  4  illustrates  how  CHAIN  works  in  practice.  Section  5  describes 
the  options  available  for  restarting  the  simulation. _ .. 
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Introduction 

A  recent  paper  (Fishman  1981)  describes  how  one  can  use  rotation 
sampling,  a  special  case  of  the  antithetic  variate  method,  to  Induce 
substantial  variance  reduction  In  the  simulation  of  a  finite  state 
Markov  chain.  Since  many  discrete  event  simulations  have  an  underlying 
Markov  structure  or  one  close  to  being  Markov,  this  variance  reducing 
proposal  has  clear  appeal.  Moreover,  for  large  and  possible  ill- 
conditioned  transition  matrices,  one  may  prefer  the  Monte  Carlo  or 
simulation  method  with  appropriate  variance  reducing  plans  to  numerical 
analysis  when  solving  for  steady-state  and  first  passage  time  distribu¬ 
tions.  In  fact,  it  may  be  the  only  feasible  method  for  some  problems. 

This  earlier  work  derives  its  cost-saving  potential  from  viewing 
the  simulation  of  k  tours  in  series  of  a  finite  (n+1)  state  posi¬ 
tive  recurrent  aperiodic  Markov  chain  as  equivalent  to  the  simulation 
of  k  replications  of  the  Markov  chain  in  parallel.  Although  the 
marginal  distributions  that  arise  with  the  two  alternative  formulations 
are  necessarily  the  same  for  corresponding  variables,  the  parallel 
formulation  all  lows  one  to  induce  joint  distributions  across  replications 
that  lead  to  a  significant  cost  saving.  The  induced  joint  distributions 
follow  from  the  use  of  rotation  sampling,  as  described  in  detail  in 
Fishman  and  Huang  (1980).  The  cost  saving  arises  In  two  ways.  Firstly, 
for  fixed  n  ,  run  time  in  the  correlated  case  Is  0(ln  k)  in  contrast 
to  0(k)  for  the  serial  simulation.  Secondly,  for  fixed  n  the 
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variance  of  an  estimator  of  interest  has  an  upper  bound  0( (1 n  k/k ) 2 ) 
for  the  correlated  case  compared  to  0(1 /k)  for  the  serial  case. 

The  present  paper  describes  an  algorithm  for  implementing  the 
essential  steps  for  k  parallel  simulations  with  correlation  along 
individual  sample  paths  induced  by  rotation  sampling.  The  paper  also 
describes  a  FORTRAN  program,  called  CHAIN,  that  one  can  use  to  perform 
the  simulation.  In  particular,  CHAIN  runs  I  independent  macroreplications 
each  of  K  correlated  parallel  microreplications  and  computes  point  esti¬ 
mates  of  interest  and  sample  variances  of  these  point  estimates. 

Section  1  introduces  the  Markov  chain  rotation  and  describes  the 
results  on  variance  reduction  derived  in  Fishman  (1981a)  for  finite  state 
chains.  For  completeness  it  also  presents  results  in  Fishman  (1981b)  for 
infinite  state  chains.  Section  2  describes  several  potential  uses  of 
CHAIN.  These  include  first  passage  time  distributions  and  steady-state 
probabilities  for  semi-Markov  processes.  Section  3  presents  Procedure  MC 
which  contains  the  essential  steps  in  carrying  out  parallel  simulation 
based  on  rotation  sampling.  It  also  describes  the  FORTRAN  CHAIN  sub¬ 
program  in  detail.  Section  4  describes  an  example  of  how  CHAIN  can  be 
used  in  practice.  The  example  is  of  the  discrete  time  Markov  chain  that 
corresponds  to  the  M/M/1  queueing  problem  with  finite  capacity  n  . 

1.  Definitions  and  Previous  Results 

Consider  a  positive  recurrent  aperiodic  Markov  chain  with  state 
space  S  *  (0,1, 2,..., n)  and  transition  probabilities  {p^;  i  ,ja0,l ,. . . ,n) 
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where  there  exists  a  positive  integer  £  <  (n-l)/2  such  that 


for  j  i-j  |  >6 


pij =  0 


i+s 

I  p.  .  =  1 

j=max(0,i-6)  1J 


It  is  convenient  to  describe  an  alternative,  but  equivalent,  representation 
to  (1)  whose  value  is  apparent  when  actually  generating  sample  paths  by 
simulation  on  a  computer.  Let  s.  denote  the  total  number  of  states  that 

J 

have  positive  transition  probabilities  from  state  j  and  let 

{m.  ;  r  *  l,...,s.)  denote  the  ordered  sequence  (m.  <  m.  r  -  l,...,s.  - 
J  •  J  J  ■  J  *  *  T  *  J 

of  the  s.  states  to  which  entry  can  occur  from  state  j  .  Then  one  has  the 

J 

representation 
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Suppose  one  wants  to  use  simulation  to  study  the  behavior  of  the 

chain  during  a  time  period  that  begins  with  exit  from  state  a  and 

ends  with  entry  to  state  b  €  S.  Consider  k  replications  of  the  simu- 

$ 

lation  experiment  run  in  parallel.  Let  denote  the  number  of 

replications  that  move  from  i  to  j  on  transition  £  and  let  K. 

J  ~ 

be  the  number  of  replications  in  state  j  at  the  end  of  transition  £.+  Let 


U,Uj ,Ug'  be  from  (1(0,1)  where  >  0  is  given.  Then  for 
J  £  , 

parallel  replications  one  can  represent  K.  .  ,  as 

JVW1 


t 

K . 


j  £ 

'  m-1  ’[Oj.r-V  V1  ^ 


(2) 


where 


and 


qj0 


=  0 


W- 1  y  •  •  •  yS  •  |  320yl  y •  •  • 
J 


U  <  X  <  V 

otherwise. 


To  ensure  that  each  replication  begins  with  an  exit  from  state  a 
and  ends  with  the  first  entry  into  state  b  ,  one  replaces 


the  original  (p^;  j=l 


’sb}  b*  Pbb  *  1  and  sets  pbfm 

br 


0  for 


The  prime  superscript  Is  used  here  for  consistency  with  the  notation  In 
Fishman  (1981a,  1981b).  However,  without  loss  of  generality,  we  take  the 
Initial  state  here  as  a  and  the  final  state  as  b  whereas  In  the  earlier 
work  the  Initial  and  final  states  were  0  and  a  respectively.  This  relabel¬ 
ing  of  Initial  and  final  states  makes  the  generality  of  the  approach  more 
apparent  to  the  reader. 
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r  *  T,...,sb  and  r  *  b  .  If  b  =  a  ,  these  modifications  are  made  after 
the  first  transition. 

If  one  uses  the  rotation  sampling  plan 


U  ■  U  + 

m  kT 

jt 


U  +  -  1 

Kjl 


0  s  U  <  1  - 


m-1 

IP 


1  -  SLjl  s  u  <  1  , 


then  the  results  in  Fishman  (1981a, b)  apply.  In  particular. 


var  0(1) 


(3) 


(4) 


and  for  ■  S  -  b  ,  the  sample  number  of  transitions  from  i  to  j 
is 


14k 


and  the  sample  mean  reward  Is 

_  1 


\  *  |j  I  1  1 

K  K  1*1  icSb  je$  U 


Then  one  has 


and,  if  jA^I  s  0(1)  , 


var  R^k  s  0( (In  k)2) 

k2  var  s  0(62(ln  k)4) 
k2  var  Rk  s  0((n6ln  k)2) 


as  n 
n  <  « 


(5a) 


(5b) 


(5c) 

(5d) 
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Here  we  speak  of  (A..)  as  the  reward  function.  Finally,  the  number  of 

•  J 

steps  to  total  absorption 


has 


and 


Tk  =  nlin(t:  Kbt  =  k) 


E  Tk  *  0(ln  k) 


var  Tk  s  0((ln  k)2) 


(6) 


Note  that  the  sampling  scheme  (3),  when  used  in  (2),  preserves  the 
identical  and  correct  probability  laws  along  the  sample  paths  of  each  of 


the  k  replications.  An  equivalent  expression  for  K^m  , 


leads  to  a  computational  saving,  is: 


which 


jr 


bry.*!  '♦ICF,5>V)  P  iQ 


mo  i-  ip  j-  iw.pi'V*)  p 


(7) 


P  =  q4.r-T 


5  *  V  V 


where 


and  x  =  x  -  LxJ  .  Here  the  cost  of  sampling  the  transitions  from  state 

j  based  on  (7),  and  using  the  inverse  transform  method*  is 

0(5.)  s  0(2*  +  1)  and  independent  of  k  whereas  sampling  cost  based  on 

(2)  is  at  best  0(k)  . 

/ 

Let  denote  the  cost  of  simulating  (1)  using  (3)  in  (7).  For  a 
countably  Infinite  state  space,  no  more  than  (26  +  1  )i  transient  states 
are  occupied  prior  to  transition  £+1  and,  therefore,  no  more  than 
(26  +  1)U2  +  l)/2  sampling  events  occur  through  transition  1  .  Since 
each  (i,j)  transition  has  cost  0(26  +1)  .one  has  s  0((6T')2)  so  that 

E  s'  £  0((6  In  k)2)  .  (8a) 

For  the  finite  state  case 

F;  S^  s  0(n  6  In  k)  .  (8b) 

These  results  compare  favorably  with  the  case  of  Independent  repli¬ 
cations  taken  In  series  where  and  R^  ,  the  sample  number  of  transi¬ 

tions  from  1  to  j  and  the  sample  mean  reward,  respectively,  have 
var  Rj^  »  0(k)  and  k  var  R^  *0(1)  .  Moreover,  for  simulation  cost  , 

one  has  0(k)  s  E  Sk  s  0(6k)  .  The  lower  bound  arises  when  n  is  small 
enough  to  allow  storage  of  all  distributions  and  their  aliases  required  by 
the  alias  method  (Walker  1977}  to  determine  the  entered  state  from  each 
existing  state  at  each  transition.  The  upper  bound  arises  from  use  of  the 
Inverse  transform  method  to  determine  the  paths.  See  Fishman  (1981b)  for 
a  more  detailed  discussion  of  this  cost. 
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It  is  also  noteworthy  that  the  desirability  of  k  parallel  replica¬ 
tions  with  rotation  sampling  relative  to  a  simulation  of  k  independent 
replications  in  series  continues  to  hold  when  {A. .}  is  not  bounded,  pro- 

*  w 

videdthat  *»>■(*„  Kjjiltl|ic5t>  <-  . 

2.  Potential  Uses 

This  section  describes  three  uses  to  which  the  previously  described 
rotation  sampling  plan  can  be  put  with  regard  to  estimation. 

First  Passage  Times 

Let 

h^  ■  probability  of  moving  from  state  a  to  state  b  for  the 
first  time  in  exactly  2.  steps 

and  (9) 

£ 

H  =  l  h.  =  probability  of  moving  from  state  a  to  state 
i=l 

b  for  the  first  time  in  no  more  than  l  steps. 


As  estimators  of 


and  one  has,  respectively. 


=  i  l 


k  ^Kb£  "  Kb,Jl-l^ 


1  r  » 

{  !  2  K, 

K  i=l  jeSb  J 


JL  =  1,2,.., 


Let  denote  the  cost  of  simulating  (1)  up  to  and  including  step 

using  the  rotation  sampling  plan  (3)  with  (7).  Theorem  1  gives  several 

*  *  , 

relevant  properties  for  h^  ,  and  Sk)l  . 


Theorem  1 .  For  a  simulation  of  k  parallel  replications  of  (1)  using 


(3)  and  (7)  one  has: 


A 


(i) 

E\  *\  • 

(ii) 

E  H  *  H  . 

?.  t 

A 

(iii) 

var  h^  <  0(  ( tS/k) 2 )  . 

(iv) 

var  0((6P./k)>)  . 

(v) 

E  SH  <  0 (  (  ;S  f  ) 2 )  . 

Proof.  Since  the  correct  probability  law  continues  to  prevail  on  sample 
paths  for  each  replication,  one  has  E  Kj^  *  k  h^  so  that  (i)  and  (ii) 

follow  immediately.  Since  var  *  0(1)  and  no  more  than  26  states 
can  transit  to  state  b  ,  it  is  clear  that 


and  that  (iii)  holds.  Since  there  are  exactly  i  steps  to  consider, 

A 

(iv)  follows  immediately.  Here  is  a  special  case  of  (5b)  with 

A..  =  1  icS.  and  A...  =  0  otherwise.  Since  exactly  f  steps  occur 

lb  b  i  J 

no  more  than  (26  +  1)U2  +  0/2  trials  are  necessary  each  with  cost  0(6) 
Therefore,  (v)  obtains. 

If  one  were  to  perform  k  independent  replications  in  series,  then 

A  A 

var  h^  =  0(l/k)  ,  var  =  OU2/k)  and  E  t  0(k)  ,  emphasizing  the 
benefit  of  rotation  sampling  plan  (3)  and  the  conciseness  of  (7). 
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Steady-State  Distribution 


p.  =  probability  of  being  in  state  j  at  an  arbitrarily  selected 

J 

time  j  =  0,1 .... ,n  . 


As  an  estimator  of  p.  one  has 


Pj  =  Gjk/Gk 


where 


*  n 


Gjk  "  k  ^t=1  ^i=0Kijt 


j=OJ  * • • » »n 


Gk*  I  Sjk 

K  j=0  JK 


Since  p.  is  a  ratio  estimator  it  has  bias 

J 


E(p  -  p.)  i  p!  'JLh. .  .co<(6ji <»  V 

J  J  J  —2  _  _  _ 

E  Gk  E  GJk  •  E  Gk 


var  G.  2  cov(G4.,G1/)  var  G. 


2  uk  c  v  jk*  k'  •  __1jk 


var  p.  *  pi  — - “ - +  - - “■ 

J  J  E  «k  ES‘ESk  E  Gjk 


and  variance 
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From  the  results  of  Section  1,  one  has 

var  s  0((6  In  k/k) 2 )  and  var  Gk  s  0((n  In  k/k) 2 ) 

so  that  both  bias  and  variance  benefit  from  rotation  sampling. 

Steady-State  Probabilities  for  a  Semi-Markov  Process 

The  results  presented  so  far  relate  directly  to  a  Markov  chain. 
However,  their  extension  to  semi-Markov  processes  is  relatively  direct 
for  the  steady-state  probabilities.  Let  y . .  denote  the  mean  time  spent 

■  J 

in  state  i  prior  to  transiting  to  state  j  .  For  the  steady-state  prob 
abilities,  one  now  has  the  estimators 

Pj  *  G^/^  j- 0,...,n 

where 

^jk  ^1=0 


and 


Then  p^  has  approximate  bias  and  variance  as  in  (13)  and  (14),  respectively, 
with  G*  replacing  Gjk  and  G^  replacing  Gk  . 
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3.  Implementation 

Procedure  MC  decribes  an  algorithm  that  one  can  use  to  encode  the 
essential  steps  for  simulating  k  replications  in  parallel  with  ini¬ 
tial  state  a  and  absorbing  state  b  .  In  practice,  one  will  want 
to  embellish  this  procedure  to  allow  for  multiple  reward  functions  and 
multiple  independent  macroreplications.  The  latter  enable  one  to  esti¬ 
mate  variances  and  covariances  of  interest. 

Figure  1  lists  a  FORTRAN  subprogram  called  CHAIN  that  enables  one 
to  simulate  an  (N  +  1)  state  Markov  chain  with  state  space  on  the  inte¬ 
gers  0,1, ...,N  using  the  parallel  rotation  sampling  plan  (7)  for  the 

purpose  of  computing  sample  mean  rewards  RPRIME(l) . RPRIME(L)  for  L 

reward  matrices  stored  in  array  A(*)  .  Each  of  I  independent  macrorepli - 
cations  begins  with  a  departure  of  K  correlated  microreplications  from 
state  INITAL  and  ends  with  the  absorption  of  all  K  microreplications  in 
state  ABSORB.  The  purpose  of  the  macroreplications  is  to  facilitate  the 
estimation  of  the  covariance  matrix  of  the  L  sample  mean  reward  functions. 
Also,  CHAIN  estimates  the  first  passage  time  distribution  from  INITAL  to 
ABSORB. 


Insert  Fig.  1  about  here. 


Procedure  MC 


Given:  a,  h,  k,  n,  (s^;  1=0 . n}  ,  {m-r:  r=1 , . . .  ,s.. ;  1=0,. ...n}  , 

{p.  :  r=l,...»si;  1=0,. ..,n}  and  (A^m  ;  r=l . s^.1  s  0,...,r 

i  »ro  .*  r  1  1  r 


1  ♦  a. 
r'  «•  0. 
K*  «■  k. 

d 


For  i=0,...,n 
For  r=1 ,. . • ,s^ 

qi,mir  " 


Initialize. 


^£=1  Pi,mi£  ‘ 


0  .  J 

5.  Go  to  8. 

6.  1  0  . 

7.  If  i=b  or  =0  go  to  23.  (Skip if  state  Is  absorbing  or  empty.) 

9.  r  «-  1.  v 

9.  Sample  U  from  U(0,1).  I 


Rotation  Sampling 


10. 

P*  -  0. 

11. 

P  -e  0. 

12. 

Q*Ki  qi.V 

13. 

q*^lqj- 

14. 

$  ■*-  Q  -  Q*  • 

75. 

X  «-  Q*  -  P*  + 

16. 

K*  *  K*  +  1 
mir  "’ir 

17. 

R'  «■  R'  +  X  A 

18. 

P*  «■  Q*  . 

19. 

F «-  $. 

(X  microreplications  move  from  1  to  m^r 
(Compute  reward.) 
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Figure  2  contains  a  partitioned  list  of  the  input  to  CHAIN.  Figure 
2a  displays  all  arrays  and  variables  for  which  an  initial  numerical  assign¬ 
ment  Is  necessary  at  the  time  at  which  CHAIN  Is  called.  Figure  2b  con¬ 
tains  all  interim  working  arrays  and  Fig.  2c  lists  all  output  arrays. 

Insert  Fig.  2  about  here. 

Note  that  the  arrays  A,  M  and  P  are  one-dimensional  in  the  subprogram 
whereas  their  counterparts  {A^}  ,  and  (P^j)  are  doubly  subscripted 

in  Section  1.  This  reduction  leads  to  a  considerable  space  saving,  especially 
when  {A..}  {m..}  and  {p..}  are  sparse  and  N  is  large.  In  terms  of 

I J  •  J  I J 

storage  space,  CHAIN  requires  0(4(ALL(L+4)  +  2L(L+3)  +  4NP  +  SIZE  +  8TT)) 
bytes  for  the  arrays  listed  in  Fig.  2.  Also,  CHAIN  has  an  upper  bound 
on  mean  execution  time  of  0(ALL  x  I  x  L  *  In  K)  . 

CHAIN  uses  the  GGUBS  random  number  generator  in  IMSL  (1977)  with  SEED 
as  the  seed  or  initial  random  number  and  returns  SIZE  uniform  deviates  on 
each  call  of  the  generator.  By  choosing  the  blocking  factor  SIZE  to  be 
large  one  reduces  the  frequency  of  calling  GGUBS,  thus  reducing  CPU  time. 
However,  the  space  requirement  for  the  uniform  deviates  is  4*SIZE  bytes. 

On  a  computer  with  limited  space,  one  may  select  a  small  SIZE  to  accomodate 
the  space  constraint.  More  generally,  an  alternative  random  number  genera¬ 
tor  can  be  substituted  by  GGUBS  with  little  effort. 

The  input  OLD  also  calls  for  explanation.  After  running  CHAIN  for, 
say  II  Independent  macroreplications  each  of  K  correlated  mlcrorepllca- 
tlons,  a  user  may  find  that  the  accuracy  of  the  sample  mean  rewards  Is 
too  low  for  Intended  purposes.  By  setting  OLD  *  II  on  a  second  run, 
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choosing  a  new  I  >  OLD  ,  using  the  original  K  and  restoring  XBAR(*) 
and  CO V(*,*)  in  the  driver  program ,  one  can  merge  the  sample  output 
from  the  first  OLD  macroreplications  with  that  from  the  subsequent 
I  -  OLD  new  macroreplications  to  produce  a  seminary  tableau.  When  this  is 
done,  the  resulting  sample  first  passage  time  distribution  is  based  on 
the  last  I  -  OLD  macroreplications  only. 

Macroreplication  versus  Microreplication 

As  Section  1  shows,  rotation  sampling  applied  to  K  parallel  micro¬ 
replications  produces  a  covariance  matrix  whose  convergence  rate  has  an 
upper  bound  0((ln  K)2/K2)  on  a  single  macroreplication.  Since  a  method 
is  not  yet  available  for  estimating  this  covariance  matrix  from  a  single 
macroreplication,  CHAIN  resorts  to  running  I  independent  macroreplica¬ 
tions  for  the  purpose  of  estimating  the  covariance  matrix  of  the  L 
sample  reward  functions.  Therefore,  the  covariance  matrix  has  a  con¬ 
vergence  rate  bounded  by  0((ln  K)2  /IK2)  and  a  run  time  of  0(1  In  K) 
for  fixed  ALL  and  L  .  Clearly  one  wants  a  K  substantially  larger  than 
I  .  In  the  example  to  be  described  next  K  *  217  and  I  ■  2 3  ,  but  other 
compromises  are  equally  reasonable. 

4.  An  Example 


To  illustrate  how  the  CHAIN  subprogram  works  in  practice,  consider 
the  Markov  chain  associated  with  the  M/M/1  queueing  problem  with  finite 
capacity  n  .  In  particular,  the  state  of  the  chain  denotes  the  number 
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of  customers  In  the  system.  In  this  queueing  problem  interarrival  times 
are  1.1. d.  from  the  exponential  distribution  with  rate  X,  service  times 
are  i.i.d.  from  the  exponential  distribution  with  rate  w  >  X  and  there 
is  a  single  server.  For  the  corresponding  Markov  chain  these  specifica¬ 
tions  imply 


poi 

=  1 

pi,i-l 

_  U) 

X+u) 

i=1 

pi,i+l 

=  _L_ 

X+u) 

i=l 

.  _x_ 

pn,n 

X+d> 

with  all  other  transition  probabilities  being  zero. 

As  objectives  consider  the  estimation  of 

W-j  =  mean  number  in  system. 

W2  =  probability  of  one  customer  in  the  system. 

W3  =  probability  of  two  customers  in  the  system. 

W4  =  first  passage  time  probability  mass  function  for  the  Markov 
chain  from  the  empty  and  Idle  state  back  to  that  state. 

Wg  =  distribution  function  associated  with  . 

Let  (Aij(£)}  denote  reward  matrix  l  and  R^(£)  denote  sample  mean 
reward  function  fc  based  on  using  (A.  .(£)}  in  (5b).  Then  one  estimates 

*  J 

Wr  W2  and  W3  by 


"2  ~  J// ''(c'  ' 

W3  -  Rj(4)/R'(2) 

<here 


Ai,1+1^  “  IS 

1*0 

1 

II 

>*|-* 

1*1 

A  0)  =  JL_ 
n  ,nv  '  A+u) 

Aio<2>  *  x 


A1,1+l^  =  A+w 

1=0 

Ai,i-1^  *  A+u 

i  =2 

An,n^  =  A*u> 

A10<«  ■  * 

A12<3»  *  TTS 

V«>  -  * 

A23^4)  *  xSi 

For  W.  and  W  we  use  the  estimators  In  (10).  These  are  computed 
**  5 

automatically  by  CHAIN. 

Figure  3  shows  a  driver  program  designed  to  Initialize  all  relevant 
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parameters  and  call  CHAIN  for  this  problem.  Here  LAM  and  W 
correspond  to  X  and  uj 


Insert  Fig.  3  about  here. 


respectively.  As  input  we  set 

LAM  =  0.9 
W  =  1.0 
N  -  49 
INITAL  *  0 
ABSORB  =  0 

I  =  23  =  8 
K  =  217  =  131072 
L  *  4 

SEED  =  1234567 
SIZE  *  40000 
TT  =  50  . 

Figure  4  shows  the  output  of  CHAIN  for  this  problem.  The  ratio 
estimators  and  together  with  their  biases  and  variances  can 


Insert  Fig.  4  about  here. 
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be  computed  by  hand  or  by  a  subsequent  subroutine  that  uses  XBAR  and 
COV  as  input  together  with  the  formulae 


E(I  .  EX,  .  EY 
EX'  EX 


[*8 


FT 


cov(X,Y)-| 
EX-EY  J 


and 


i  E2Y  rvar  X  2  cov(X.Y)  .  var  Y  , 
var(Y/X)  -  pY  '-~£rx - F7  .  fv  +  ~T7~T  J 


EX  •  EY 


We  have  chosen  to  do  them  by  hand.  They  are 


W1  =  8.7396 


W2  =  .0905 


E  Y 


W3  *  .0814 


E(W1-W1)  =  -.1037  x  10-4  E(W2-W2)  =  .4304  x  10"  E(W3-W3)  =  .3735  x  10' 
var  W,  =  .2646  x  10"3  var  W2  =  .3780  x  10"8  var  W3  =  .2936  x  10' 


From  Gross  and  Harris  (1974,  Sec.  2.5)  one  has 

W1  =  8.7410  W2  =  .0905  W3  =  .0814  , 

which  confirm  the  performance  of  CHAIN. 

All  computing  was  performed  on  an  IBM  370/155  computer.  For 
FORTRAN  G  (level  21),  the  CHAIN  Program  requires  8592  bytes  of  space. 
For  FORTRAN  H  (level  21.8)  it  requires  7380  bytes. 
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5.  Restarting  the  Simulation 

Situations  will  occur  in  which  a  user  of  CHAIN  does  not  have  a  good 
a  priori  estimate  of  the  running  time  or  the  statistical  reliability 
to  be  expected  from  a  specified  set  of  1^  macroreplications  each  of 
microreplications.  At  least  two  alternatives  exist  for  dealing  with  this 
case.  A  user  may  make  a  preliminary  run  and  determine  that  I2  additional 
macroreplications,  each  of  K2  =  microreplications,  are  necessary  to 
achieve  the  desired  accuracy  within  budget.  CHAIN  is  designed  to  accommodate 
this  alternative  and,  provided  that  the  sample  means  and  covariances  accumu¬ 
lated  on  the  first  1^  macroreplications  are  restored  prior  to  collecting  the 
additional  I2  macroreplications,  the  routine  prints  the  global  sample 
means  and  covariances  at  the  end  of  all  1^  +  I2  macroreplications. 

The  second  alternative  arises  when  on  the  basis  of  the  first  run  a 
user  determines  that  I2  macroreplications  each  of  K2  *  K1  microrepli¬ 
cations  is  the  most  desirable  way  of  achieving  the  desired  accuracy.  Here 
the  restoration  feature  of  CHAIN  does  not  apply.  Let  J  and  i  be  the 
sample  mean  vector  and  sample  covariance  matrix  of  Y  obtained  on  the 
first  run  with  1^  and  K-j  .  Let  Y  and  n  denote  the  corresponding 
quantities  on  the  second  run  with  I2  and  K2  .  Then  the  overall  sample 
mean  vector  is 

I '  TITT-IiT  (11  I  +  T> 

and  its  sample  covariance  matrix  is 
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A1 though  the  computation  of  I  and  r  are  not  features  of  CHAIN 
one  can  easily  add  them  if  desired.  However,  as  in  the  case  of  the 
ratio  estimators  in  Section  4,  the  relative  importance  of  this  feature 
will  differ  from  user  to  user. 
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SUBROUTINE  CuAi  N  (K.OLD,  I,  NP.INITAL,  A BSORB ,S, H.SUHS, ALL,  P. 

1  KPRIM£,KSTAB,Q,L,ASIZE,A,RPRIHI, XBARaCOV,COEFP, 

2  TT,FRBAB,FRVAB,CUHBAR,CUHVAR, SEED, SIZE, V) 

C  •••  RKOGRAH  FOB  SIMULATION  OF  MARKOV  CHAIN  IXTH  N+1  STATES 

C  USING  DOTATION  SAMPLING. 

INTEGER  HP,ABSOBB,  ALL,  B,  COO  NT,  I  ( IN  IT  AL,  ISEED,  IT,J,KPRIHE{NF), 

1  KSTAR  (HP)  ,L,LA,LI,LIJ,K,LL,  LM,R, H  (ALL)  ,  N,  OLD,  OR  EGA,  Z, 

2  PSTAB, QSTA8. OMEGA, SEED, SIZE, S(NP) ,S J. S A , SS , SOBS (NP) ,TT, 

3  ASIZE 

BE1L*4  A  (ASIZE)  , TRCV ,  CUHCV, V  (SIZE) 

REAL'S  C,  CC ,  COEP F  (L)  ,  CO V  (L,  I)  >  P (  ALL)  ,  PP#  PB  AB«  Q  (ALL) ,  QQ,  QB  AB, 

1  APBIHE(L)  ,RSE£D,XBAR  (1)  .  0,PREQ,COi.  PBVAB  (TT)  ,  PBBAI  (TT|  , 

2  CUHbAB  (TT)  , CUH VAR  (TT)  . FP, FPV AR , REG 


DESCRIPTION  OF  VARIABLES  - 


•  •• 

•  •• 

•  ••  l  (•)  s  REWARD  MATRIX 

•••  ABSORB  *  ABSORBING  STATE 

•••  ALL  =  SOB  OF  S(J)  FOR  ALL  J 

ASIZE  =  SIZE  OF  A  ARRAY 

•••  B  *  TEST  VARIABLE 

•••  C  «  ((L1-1)/LI) 

•  ••  CC  *  l  (LI-OLD-1)  /  (LI-OLD)  ) 

*••  COEFF (•)  *  COEFFICIENTS  OF  VARIATION 

COUNT  *  UNIFORH  DEVIATE  COUNTER 
•••  CO V (• , •)  *  COVARIANCE  HATRII 

CUH  =  ACCUMULATION  OF  FREQ 

•••  CUHDAR (*) *  MEAN  OF  CUH 
C  •••  cuncv  =  COEFFICIENT  OF  VARIATION  FOR  CUH 
C  •••  CUHVAM*)*  VARIANCE  OF  CUH 
c  •••  FF  *  SAflPLE  BEAN  FOR  FIRST  PASSAGE 

C  •••  FFVAR  *  SAHPLE  VARIANCE  FOB  FIRST  PASSAGE 

C  •••  FBBAB (*)  *  MEAN  OF  FREQ 

C  •••  FHCV  *  COEFFICIENT  OF  VARIATION  FOR  FREQ 

C  •••  FREQ  *  NUMBER  OF  NEk  TRANSITIONS  INTO  ABSOIRINC  STATE 

C  •••  FRVAR (•)  «  VARIANCE  OF  FREQ 

C  •••  I  *  DESIRED  RUBBER  OF  HACROREPLICATIONS 

C  •••  INITAL  *  INITIAL  STATE 

C  ISEED  =  INITIAL  SEED 

c  •••  IT  X  TRANSITION  COUNTER 

C  •••  J  *  INDEX  FOR  STATES 

L  •  K  *  NUMBER  OF  PARALLEL  HICROIEPIICATIONS 

C  •••  KPBI HE  (•) *  NUMBER  IN  STATE  AT  END  OP  IIANSITION 

C  •••  KSTAR  (*)  *  NEXT  KPRIHE  (*) 

C  •••  L  *  NUMBER  OF  DESCRIPTORS  TO  BE  ESTIHATED 

C  LA  *  LL  4  1 

C  •••  LI  •  INDEX  FOR  I 

■  •••  LI J  *  1  +  OLD 

•••  LJ  »  INDEX  FOB  STATES 


00000100 
00000110 
00000120 
00000130 
000001*0 
00000  ISO 
00000160 
00000170 
00000180 
00000190 
00000200 
00000210 
00000220 
00000230 
00000240 
000002SO 
00000260 
00000270 
00000280 
00000290 
00000300 
00000310 
00000320 
00000330 
000003*0 
000003S0 
00000360 
00000370 
00000380 
00000390 
00000400 
00000410 
00000420 
000004 30 
00000440 
00000*50 
00000*60 
00000*70 
00000*80 
00000*90 

ooooosoo 

00000510 

00000520 

00000530 

000005*0 

00000550 

00000560 

00000570 

00000580 

00000590 


Fig.  1.  CHAIN  Subroutine 


r n  r.  n  r  n  r.  r,  n  r.  r.  r,  r>  r>  n  r.  n  r.  r.  n  n  n  n  r,  n  n  n  r.  n  r 
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*•*  11  »  INDEX  FOE  L 

•••  LB  *  INDEX  FOR  L 

•••  H{*|  =  STATES  THAT  CAM  3E  ENTERED  PBOM  J- 1 

•••  M  *  NUMBER  OF  HIGHEST  STATE 

•  ••  Mi  *  BUNDER  Of  STATES  (N+ 1| 

•••  OLD  »  NUNBE8  OF  B ACROR EPL1CAT1QNS  ALBEADT  PEI FORHEO  (OLD<I) 

•••  OMEGA  «  TEST  VARIABLE 

•••  P(»)  »  TB ANS1T10N  BATHII 

•  ••  tp  *  uo  PBOB  TIB E  BEPOBE 

•••  PBAB  *  PBACT10BA1  PART  OP  PP 

•••  PSTAB  »  INTEGER  PART  OP  PP 

Q(*)  *  P BOH  P 

•••  UQ  *  FROM  U  AMD  K 

JBAR  *  F b  ACT  IC N AL  PAST  OP  QQ 

•  ••  JSTAR  *  INTEGER  PART  OP  30 

•••  B  «  INDEX  FOB  STATES 

•••  REG  *  ACCUMULATED  FREG  POE  TRANSITIONS  CREATES  THAR  TT-1 

•  ••  BPB I  HE  ( •)  *  CUMULATIVE  BENABD  POR  ESTIMATORS 

•••  RSEED  *  RANDOM  CENLBATOR  SEED,  SEAL  VALUED 

•  SEED  «  RANDOM  GENE  BATOR  SEED 

•••  SIZE  «  bLOCK  SIZE  POR  RAMDOn  RUMBER  GENERATOR 

•••  S(J)  «  NUH3LB  OP  STATES  THAT  CAB  BE  BVTEIED  PBOH  J-1 

•••  SJ  =  S(J) 

•••  SB  «  SUHS(J) 

•••  SS  *  MIN  { NUMBER  OP  TRANSITIONS ,  TTJ 

•  ••  SUMS  (J)  =  SUM  OF  S  (1)  1 HBU  S  (J- 1)  ,  AMD  SUBS|1)>0 

•••  TT  *  MAXIMUM  NUMBER  OP  TRAMSITXOBS  TO  LOOK  AT 

•••  0  *  CURRENT  UNIFORM  DEVIATE 

•••  V(»J  -  UNIFORM  DEVIATE  ARBAT 

•••  X  <  NUMBER  OP  TRANSITIONS  IM  CUBBBBT  MOVE 

•••  XBAB  (•)  «  MEAN  MATRIX 

C  •••  INITIALIZE  VARIABLES. 

PP  »  0.000 
PPVAR  *  O.ODO 
SS  =0 
ISEED  »  SEED 
RSEED  *  SEED 
COUNT  =  SIZE 
DO  80  J=1,TT 

PRVAR(J)  *  O.ODO 
FB3AR  (0)  *  O.ODO 

CUnVAR  (J)  a  O.ODO 
80  COKBAR  (J)  ->  O.ODO 

DO  100  J«1,NP 

SK  a  SOMS(J) 

Ct<SK+l)  «  P(SK+1) 

IF  (S(J).L1.2)  GO  TO  100 
SJ  as  (J) 


00000600 
00000610 
00000620 
00000630 
000006*0 
000006 SO 
00000660 
00000670 
00000680 
00000690 
00000700 
00000710 
00000720 
00000730 
000007*0 
00000750 
00000760 
00000770 
00000780 
00000790 
00000800 
00000810 
00000820 
00000830 
000008*0 
00000850 
00000860 
00000870 
00000880 
09000890 
00000900 
00000910 
00000920 
00000930 
000 009*0 
00000950 
00000960 
00000970 
00000980 
00000990 
00001000 
00001010 
00001020 
00001030 
000010*0 
00001050 
00001060 
00001070 
00001080 
0 0001090 


Fig.  1  ( Continued ) 
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00  90  B*2,SJ 

00001100 

90 

0  (Sr.+  h)  «  Q(SK*R-1)fP(SKtB> 

00001110 

u  (SK  +  SJ)  »  1.000000 

00001120 

100 

CONTINUE 

00001130 

c  ••• 

CliECK  IF  X  BA  R  AND  COV  ABE  BESTOBEO. 

00001140 

IF  (OLD.GT.O)  GO  TO  115 

00001150 

00  110  11*1, L 

oooomo 

X BAB  (Ll|  -  0.000 

00001170 

COEFF(LL)  »  0.000 

000011B0 

DO  110  Lfl-  1.L 

00001190 

1  10 

COV  (LL.Lfl)  *  0.000 

00001200 

1  15 

LI J  «  1  +  OLD 

00001210 

DO  1 18  LL«1,L 

00001220 

X BAR  (LL)  =  NBAS  (LL)  *OLD 

00001230 

DO  118  LB*1,L 

00001240 

116 

COV  (LL. Lfl)  *  COV (LL. Lfl) *OLD 

00001250 

00001260 

C  ••• 

START  MAIN  LOOP. 

00001270 

DO  540  L1*L1J .1 

00001280 

c  ••• 

INITIALIZE  VARIABLES  FOB  THIS  REPLICATION. 

00001290 

B  =0 

00001300 

C  =  (LI-  1.  ODO)  /LI 

00001310 

Cc  ■=■  ( LI-OLD- 1.  ODO)  /(LI-OLD) 

00001320 

CUB  =  0.  ODO 

00001330 

BEG  *  0.  ODO 

00001340 

IT  =  1 

00001350 

OB  EG  A  -  0 

00001360 

DO  no  LL-  1.L 

00001370 

no 

RPKIBE  (LL)  =  0.000 

00001380 

DO  130  J*  1 . NP 

00001390 

KSTAF/J)  *  0 

00001400 

130 

KPhIBE{J)=  0 

00001410 

C  ••• 

START  THIS  REPLICATION  HITH  ALL  BICIOBEPLICATIONS  IN 

INITIAL  STATE00001420 

J  =  INI TAL  +  1 

00001430 

KPB1HE  (J)  =  K 

00001440 

KSTAR(J)  «  K 

00001450 

c  ••• 

LOOK  AT  INITIAL  STATE. 

00001440 

J  *  1NITAL  *  1 

00001470 

GO  TO  400 

00001480 

c  ••• 

LOOK  AT  THE  NEXT  STATE. 

00001490 

300 

J  «  J+1 

00001500 

c  ••• 

SKIP  THE  ABSORBING  STATE. 

00001510 

IF  (J.  Ey.  ABSORB*  1)  J  ■  J*1 

00001520 

c  ••• 

IF  NONE  IN  THIS  STATE,  LOOK  AT  THE 

NEXT  STATE. 

00001530 

400 

IF  (KPRIBE(J) .Eg.O)  GO  TO  300 

00001540 

c  ••• 

FIND  THE  NUHBER  OP  STATES  THAT  CAR 

BE  ENTUSO  FEOfl  STATE  J-1. 

00001550 

SJ  *  S(J) 

00001560 

C  ••• 

START  UlTn  T HE  FIRST  STATE  THAT  CAN 

BE  ENTERED  FEOfl 

STATE  J-1. 

00001570 

R  *  1 

00001580 

c  ••• 

POINT  TO  THE  NEXT  DEVIATE  TO  USE. 

00001590 

Fig.  1  ( Continued ) 
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c  ••• 

c  ••• 

c  ••• 
c  ••• 

*20 

•  •• 
c  ••• 

C  ••• 

500 

c  ••• 

C  ••• 

c  ••• 

c  ••• 
c  ••• 
c  ••• 

510 

c  ••• 

c  ••• 
c  ••• 
c  ••• 
c  ••• 

c  ••• 


COUNT  *  COUNT  •»  1 

GST  NEN  ABBA*  Or  DEVIATES  IF  NEEDED. 

IF  (COUNT. LE. SIZE)  GO  TO  *20 
COUNT  *  1 

CALL  BANDOn (SEED. V, SIZE) 

CALL  GGUBS (RSEED.SIZE,  V) 

GET  THE  DEVIATE  POINTED  TO  BI  COUNT. 

TEANSrOBB  DEVIATE. 

U  «  DBLE  ( ABOD  (KPBIBE  (J)  *V  (COUNT)  ,  1. )  ) 

NAITE  (3.  20*0)  U ,  V  (COUNT) 

TBANSFEB  ALL  OUT  Or  THIS  STATE  FOB  TRANSITIONS. 

KSTAB(J)  «  XSTAB(J)  -  KPBIBE  (J) 

INITIALIZE  QQ.  QSTAR,  AND  QBAB. 

VV  >0.0 
QSTAR  *  0 
QBAB  =0.0 

SAVE  THE  LAST  OCCUBBENCE  OF  Qtf,  QSTAB,  AND  QBAB. 

pp  =  qw 

PSTAB  =  QSTAB 

PDAB  =  QSAB 

COBPUTt  Nth  VALUES  FOR  Gw,  QSTAB,  AND  QBAB. 

QQ  »  KPFIHE  (J)  *w  (Suns  (JJ  +  B) 

QSTAB  =  I  DINT  (QQ) 

QBAB  «  DPIODU'Q,  I.ODO) 

FIND  THE  «  CF  TRANSITIONS  TO  THE  B-TH  STATE  THAT  CAN  BE  ENTERED. 

FROtt  J-1. 

X  =  (wSTAR-PSTAR)+.5*(DSIGN(1.0D0,O-PBAB)+ 

1  DSIGN( I.ODO, QBAF-0) ) 

ADD  THESE  TaAKSITIONS  TO  THE  ENTERED  STATE  OCCOPANC!  VECTOR. 

KSTAR(fl(SUHS(JHN)  +  1)  *  KSTAR  (H  (SUNS  (J)  f  B)  -f  1 J  +X 
FIND  THE  NUHBER  CF  TRANSITIONS  BADE  SO  FAN  FOB  THIS  STATE. 

II  *  D-»X 

ACCUHULATE  THE  BENABDS  FOR  THIS  STATE  TtANSITION. 

DO  5.0  LL-  1 , L 

IP  (A(ALL*(IL-1)+SUHS(J)+R).EQ.O.O.OB.I.EQ.O)  GO  TO  510 
BPBIHE(LL)  *  RPHIH£(LL)+A(ALL*(LL-1)+SOBS(J)+B|*I 
CONTINUE 

GO  TO  THE  NEXT  STATE  THAT  CAN  BE  BEACHED  PBOH  J-1. 
h  =  B+  1 

IT  NOT  ALL  TRANSITIONS  HERE  HADE  FOB  THIS  STATE,  TBI  AGAIN. 

IF  (D.LT.RPRIHEfJ) )  GO  TO  500 

ACCUHULATE  THE  RUBBER  OP  TRANSITIONS  BADE  SO  FAB  FOB  ALL  STATES. 
ONEGA  »  OHEGA+B 

CLEAR  THE  NUflBEA  OF  TRANSITIONS  FOR  THE  STATE  COUNTER. 

B  =0 

IF  NOT  ALL  TRANSITIONS  HERE  BADE  VET,  TBT  AGAIN. 

IF  (OHECA.  LT.  K-KPRI  HE  (ADSORB*  1))  GO  TO  300 
ONEGA  =  0 

ACCUHULATE  FB3AR,  FNVAB,  CUHbAR  AND  CUBVAB. 

FKEw  *  KSIAR  (A  BSORB+ 1)  - KPB1  HE ( ABSOBB+ 1| 


00001600 
00001610 
00001620 
00001630 
000016*0 
00001650 
00001660 
00001670 
00001680 
00001690 
00001700 
00001710 
000017 20 
00001730 
000017*0 
00001750 
00001760 
00001770 
00001780 
00001790 
00001800 
00001810 
00001820 
00001830 
000018*0 
00001850 
00001860 
00001870 
00001880 
00001890 
00001900 
00001910 
00001920 
00001930 
000019*0 
00001950 
00001960 
00001970 
00001980 
00001990 
00002000 
00002010 
00002020 
000020  JO 
000020*0 
00002050 
00002060 
00002070 
00002080 
00002090 


Fig.  1  ( Continued ) 
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1 1  (IT.Ew.1J  FEEw  *  KSTAB  (ABSORB4  1) 

00002100 

FF  *  FI  4  IT*FBE3 

00002110 

FFVAP  *  FFVAB  4  IT**2*FBEU 

00002120 

CUM  *  CUM  4  FREQ 

00002130 

LL  •  NINO  (TT, IT) 

000021*0 

IF  (IT. LT. TT)  CO  TO  515 

00002150 

BE 0  -  RED  4  FREQ 

00002160 

IF  (KSTAB (AbSONBft) .LT.K)  00  TO  522 

00002170 

FBEQ  «  REG 

00002180 

515 

FBBAR(LL)  *  FRBAR  (LL)  4  FREQ 

00002190 

CUHBAK  (LL)  *  CUHUAB  (LL)  4  CUM 

00002200 

F HYAR(LL)  «  FRVAB (LL)  4  PBEQ**2 

00002210 

CUHVAB  (LL)  *  CUHVAB  (LL)  4  C0H*»2 

00002220 

C  ••• 

BESET  ThE  STATE  OCCUPANCY  VBCTOB. 

000022JO 

522 

DO  525  LJ>1,NP 

000022*0 

KPUIHE(LJ)  »  KSTAB  (LJ) 

00002250 

SJ  «  S  (LJ) 

00002260 

DO  525  B*  1 , S J 

00002270 

525 

KPRIRE(H  (SUMS  (LJ)  48)  41)  «  KSTAB  (B(SUHS  (LJ)  4B)  4U 

00002280 

C  ••• 

START  NEXT  TRANSITION  IN  STATE  0. 

00002290 

J  *  0 

00002300 

c  ••• 

COMPUTE  TUE  NUMBER  OF  TRANSITIONS  BADE  SO  FAB. 

00002310 

II  *  IT  4  • 

00002320 

c  ••• 

IF  NOT  ALL  ABSORBED,  TRY  AGAIN. 

00002330 

IF  (KPRIHE (ABSOHB4 1) . LT. K)  GO  TO  300 

000023*0 

c  ••• 

RECURSIVE  COMPUTATIONS. 

00002350 

IF  (LL.GE.TT)  GO  TO  529 

00002360 

c  AM 

ACCUMULATE  CUHBAR  AND  CUHVAB  FOB  •  OF  TRANSITION  STEPS  >  TT- 1. 

00002370 

LL  =  LL  4  1 

00002380 

DO  527  LH=LL,TT 

00002390 

CUNVAR(LH)  *  CUHVAB (LB)  4  C0M**2 

00002*00 

527 

CUHBAR  (LM)  *  CUHBAR(LH)  4  CUM 

00002*10 

C  ••• 

COMPUTE  THE  COVARIANCE  MATRIX  RECURSIVELY. 

00002*20 

529 

IF  (LI. Eg.  1)  GO  TO  535 

00002*30 

DO  530  ll-I.L 

00002**0 

DO  530  LM«1,l 

00002*50 

530 

COV  (LL,  Lfl)  »  (  (LI-2)  *COV  (LL,  LB)  4  (EFBIBE  (LL)  /K-XBAR  (LL)  )  • 

00002*60 

00002*70 

00002*80 

00002*90 

00002500 

00002510 

00002520 

00002530 

000025*0 

00002550 

00002560 

00002570 

00002580 

00002590 


1  (B  PBI BE  { LH)  /K-X BIB  (LB)  )  *C)  /  (LI-  1) 

C  •••  COMPUTE  ss. 

535  SS  *  HIN0(TT,HAX0(SS,XT-1)) 

C  •••  COMPUTE  TUE  SAMPLE  MEAN  VECTOR  RECURSIVELY. 

DO  5*0  tl»  1,L 

XBAR(LL)  »  C*XBAR  (LL)  4  ( 1.  ODO-C)  •RPRXHE  (LL)/K 
5*0  CONTINUE 

C  •••  END  MAIN  LOOP,  COMPLETED  1  BACBOBEPLXCATXOBS. 

C  ••• 

C  •••  COMPUTE  FBBAS,  FRVAR,  CUHUAB,  AA'D  CUHVAB  DIBECTLT. 

IP  (I.NE.1)  am  «  (K**2.0D0)  *(I**2.0DO)*(1-1.0DO) 

DO  S*5  LL-1.TT 


Ffg.  1  ( Continued ) 


IP  (Z.  Ew.  1)  CO  TO  544 

00002600 

PRVAR(LL)  -  (l»PRVAR (LL)-PRBAR(LL) ••2)/REG 

000026 tO 

CURVAR(LL)  *  (I*CURVAR(LL) -CURBAR (LL) ••2J/REG 

00002620 

544 

CURBAR  (LL)  *  CURBAR  (LL)  /  (K»I ) 

00002630 

S«5 

FRBAR(LL)  -  PRbAR  (LL)/(K*X) 

00002640 

00002650 

DO 

550  LL*1,L 

00002660 

DO  550  Lfl* 1, L 

00002670 

550 

COV  (LL.LR)  «  COV  (LL,  Lfl) /L 

00002680 

00 

570  LL*  1.L 

00002690 

IP  (XBAF  (LL) .  NE.  0)  COAPP(LL)  *  DSQBT  (COV  (LL.LL) )  /(BAR (LL) 

00002700 

IP  (Ll.EQ.  1)  GO  TO  570 

00002710 

LA  «  LLt 1 

00002720 

DO  560  Lfl* LA, L 

00002730 

IP  (COV  (LL.LL)  *COV  (LR,LH)  •  GT.  0) 

00002740 

1 

COV  (LH, LL)  *COV  (LL.LR)  /  DSQBT  (COV  (LL.LL)  «OOV  (LR.LR)  ) 

00002750 

560 

CONTINUE 

00002760 

570 

CONTINUE 

00002770 

00002780 

c  +•• 

PRINT  RESULTS. 

00002790 

00002B00 

SEED  »  RSEED 

00002B 10 

WRITE  ( 3, 1000)  NP.IN1TAL, ABSORB, ALL.L.K, I, ISEED.SEED. SIZE 

00002820 

1000 

POUR AT (1  Hi. 'RESULTS  OP  ROTATION  SARPLIN6 

FOB  A  HARKOV  CHAIN',///, 

00002830 

1* 

NO.  OF  STATES 

—  .110//. 

00002840 

2* 

INITIAL  STATE 

-.110//, 

000028S0 

3* 

ABSORBING  STATE 

-.no//. 

00002860 

4* 

TOTAL  NO.  OF  (I,J)  PAIRS 

-.no//. 

00002870 

5* 

NO.  OF  DESCRIPTORS 

-.no//. 

00002880 

6* 

NO.  OF  CORRELATED  HICROREPLICATIONS 

—.no//. 

00002890 

7* 

NO.  OP  INDEPENDENT  R ACI OR BPL ICA TJ 0N5 

— ,310//, 

00002900 

8* 

INITIAL  SEED 

••.no//. 

00002910 

9* 

FINAL  SELD 

—  .110//, 

00002920 

C' 

BLOCKING  FACTOR 

—.no/) 

00002930 

BEHIND  11 

00002940 

WRITE  (3.20001 

00002950 

WRITE  (3,2030) 

00002960 

WRITE  (3,2040)  (XB  AR  (LL)  ,  LL*  1  .  L) 

00002970 

WRITE  (3.2010) 

00002980 

WRITE  (3.2030) 

00002990 

WRITE  (3,2040) (COEPP(LL) ,LL*1,L) 

00003000 

WRITE  (3,2020) 

00003010 

WRITE  (3.2030) 

00003020 

DO 

560  LL* 1 , L 

00003030 

WRITE  (3,2045)  LL,  (COV  (LL.LR)  ,  Lfl*  1,1) 

00003040 

580 

CONTINUE 

00003050 

2000 

PORRAT (//, •  SAHPLE  BEAN  VECTOR',/, 

00003060 

1 

00003070 

2010 

PORRAT (//,'  SAHPLE  COEFFICIENTS  OP  VARIATION',/. 

00003060 

1 

►•••.,/) 

00003090 

Fig.  1  { Continued ) 


2020 

2030 

2040 
204 It 


610 
o  20 

3000 

3010 

3020 


3030 

4000 
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FGHH AT  (//,  •  SAMPLE  COV  A&I  ANCE/COIBELATIOM  HATHA',/, 

I  . * . . . 

FOBHAT(9X,« 1',14X,'2',14X,'3' , 1**,' *• , I*!.' 5* , UX, *6' , 141,' 7* , 
1  14X,'8'.14X.'9 *,//) 

FOBH  AT  (2X, 9  (D15.  8)/) 

FOB  BAT  (12,9  (015.  8) /| 

FF  «  rP/(K*I) 

FFVAR  *  FFVAB/ (K*I)  -  FF»*2 

00  620  LL* 1, SS 
FHCV  =0.0 

cuncv  =0.0 

IF  (FRBAR  (LL)  .  HE.  0.  000)  PRCV*DSQRT  (FRVAR  (LL)  )/FBBAR(LL) 

IF  (CUBBAR  ILL) . BE. 0.000)  CUHCV'DSUBT (CUBVAB (LL) ) /CURB AH (LL) 
IF  (ROD (LL, 50 )  .  ME.  1 )  GO  TO  620 
WRITE  (3,3000) 

IF  (LL.GT.  1)  GO  TO  610 
WAITE  ( J, 3010)  ff, FFVAR 
WHITE  (3,3020) 

■  CITE  (3,3030)  LL.FBBAH (LL) , FHVAB (LL) , PRC  V,  CO  MBA  B  (CL)  , 

1  CURVAK(LL) ,CUHCV 

IF  (OLD.EvJ.O)  PETUHM 


LL  =  I  -  OLD 
WRITE  (3,4000)  LL 

FOSR AT (1h 1, 'FIRST  PASSAGE  TIRE  (T)  DISTRIBUTION' ,/, 

1  *  . . . . 

FORMAT (/, lx, 'SAMPLE  MEAN (T)  =',015.8,51, 

1  'SAMPLE  VARIANCE  (T)  =’,015.8) 

FORMAT (/,33X, 'HASS  FUNCTION* , 42X, 'D2STBIBDT20N  FUNCTION',//, 

1  351, • VARIANCE* , 10X, • COEFFICIENT* , 

2  31X,' VARIANCE', 10X, 'COEFFICIENT',/, 

3  •  I',10X,*PB  (T=IJ ' ,  11X, 'OF  PR(T=I)  ',91,'Or  VARIATION'. 

4  1 0X, *  PR (T<*I) ' , 10X, *  OF  Pi(T<-I) ' ,91, 'OF  VARIATION', 

5  . . ,9X,' - ',  111,* - '.91. • - ', 

6  10X,  - -  ,  101,  •  — — - - - - •) 

FORM AT ('  ',15,6(51.015.8)) 

POSnAT (/, •  FIRST  PASSAGE  DISTRIBUTION  IS  BASED  ON  THE  LAST 
1  15, »  RACBOREPLICATIONS. •) 

RETURN 

END 


00003100 
00003110 
00003120 
00003130 
00003140 
00003150 
00003160 
00003170 
00003180 
00003190 
00003200 
00003210 
00003220 
00003230 
00003240 
00003250 
000032 60 
00003270 
00003260 
00003290 
00003300 
00003310 
00003320 
00003330 
00003340 
00003350 
00003360 
00003370 
000033B0 
00003390 
00003400 
00003410 
00003420 
00003430 
00003440 
00003450 
00003460 
0000347(1 
00003480 


Fig.  1  ( Continued ) 
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Fig.  2  Input  to  CHAIN  Subprogram 
(a)  Data  Input 


VARIABLE 

TYPE 

DESCRIPTION 

A 

REALM(ASIZE) 

A(ALL*(LL-1)  ♦  SUHS(J)  ♦  R)  ■  reward  received  when 
a  replication  junps  from  state  J  -  1  to  state 
M(SUMS(j)  ♦  R)  for  R  »  1,...,S(0)  for  reward 
function  LL  »  1,2 . L 

ABSORB 

INTEGER 

Absorbing  state 

0  s  ABSORB  s  N  ♦  1 

ALL 

INTEGER 

Total  number  of  arcs  -  SUKS(NP)  ♦  S(NP) 

ASIZE 

INTEGER 

ALL*L 

I 

INTEGER 

Desired  number  of  Independent  macro replications 

INITAL 

INTEGER 

Initial  state 

0  s  INITAL  s  N  ♦  1 

K 

INTEGER 

Number  of  parallel  mtcrorepllcatlons  per 
mdcrorepl Icatlon 

L 

INTEGER 

Total  maaber  of  reward  functions 

M 

INTEGER(N+1) 

N($UM{J)  ♦  LR)  ■  LRth  of  S(D)  states  to  which  a 
replication  can  move  from  state  J  -  1 

LR  -  1,...,S(0) 

NP 

INTEGER 

NP  ■  N  ♦  1  ■  total  number  of  states 

OLD 

INTEGER 

If  OLD  *  0  simulation  proceeds  to  run  I  macro- 
replications 

If  OLD  >  0  simulation  proceeds  to  run  I  -  OLD 
additional  macrorepll cations 

P 

REALM  (ALL) 

P(SUM(J )  *  LR)  ■  probability  of  moving  from  state 

J  -  1  to  M(SUM(J)  ♦  LR)  LR  -  1 . S(J) 

S 

1 

INTEGER(N+1 ) 

S(J)  «  number  of  sutes  that  can  be  entered  from 
state  J  -  1 

!  SEED 

INTEGER 

Initial  value  for  random  number  generator 

|  SIZE 

INTEGER 

Each  call  to  the  random  number  generator  returns  a 
block  of  SIZE  uniform  deviates 

! 


(continued) 
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Fig.  2  (Continued) 
(a) 


VARIABLE 

TYPE 

DESCRIPTION 

SUNS 

INTEGER(NP) 

SUNS(J)  -  0  J  ■  1 

J-1 

-  I  S(I)  J  -  2 . NP 

1-1 

TT 

INTEGER 

Number  of  cells  In  sample  first  passage  time  distri¬ 
bution.  Last  cell  estimates  probability  of  absorp¬ 
tion  at  time  a  TT 

[  (b)  Working  Arrays 


VARIABLE 

TYPE 

DESCRIPTION 

KPRIME 

INTEGER(NP) 

Distribution  of  mlcrorepllcatlons  by  state  at 
end  of  a  transition 

KSTAR 

1 

INTEGER (NP) 

Distribution  of  mlcrorepllcatlons  by  state  at 
beginning  of  a  transition 

Q 

REAL*8(ALL) 

Q(SUN(0)  ♦  LR)  •  probability  of  moving  from 
state  J  -  I  to  state  N(SUN(J)  ♦!), 

N($UM(J)  ♦  2),...  or  H(SUH(J)  ♦  LR) 

LR  •  1,...,S(J) 

RPRIME 

REAL*8(L) 

Accumulated  rewards  for  L  reward  functions 

V 

REAL*4(SIZE) 

Space  to  store  uniform  deviates 

(continued) 
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Fig.  2  (Continued) 


(c)  Arrays  Used  to  Summarize  Data  on  I  Macro  replications 


VARIABLE 

— 

TYPE 

DESCRIPTION 

COEFF 

REAL*8(L) 

COEFF(JI)  ■  sample  coefficient  of  variation  of 
XBAR(JI)  01-1 . L 

COV 

R£AL*8(L,L) 

At  completion  COV (01.  02)  contains  the  sample 
covariance  of  XBAR(01)  and  XBAR(02)  for 

J2-01 . L  and  Jl-1 . L  and  C0V(J2,J1)  con¬ 

tains  the  sample  correlation  between  XBAR(01) 

and  XBAR(02)  for  01-1 . 02-1  and 

02-2, ...,L  . 

CUNBAR 

REAL*8(TT) 

CUMBAR(Ol)  -  sample  probability  that  absorption 
occurs  on  a  step  s  01  for  01-1 , . . . ,TT  -  1 ; 
CUMBAR(TT)  -  1 

CUMVAR 

REAL*8{TT) 

CUMVAR(Ol)  -  Sample  variance  of  CUMBAR(Ol) 

01-1 . TT 

FF 

REAL*8 

Sample  man  of  first  passage  time. 

FFVAR 

REAL*8 

Sample  variance  of  first  passage  time. 

FRBAR 

REAL*8(TT) 

FRBAR(Ol)  ■  sample  absorption  probability  at 
step  01  for  01-1,..., TT  -  1 

TT-I 

FRBAR(TT)  •  1  -  )  FRBAR(Ol) 

01-1 

FRVAR 

REAL*8(TT) 

FRVAR(Jl)  -  sample  variance  of  FRBAR(OI) 

01-1,..., TT 

XBAR 


REAL*8(l) 


XBAR(OI)  ■  the  staple  man  reward  for  reward 
function  01-1,. ..,L 


INTEGER  X. J.  UPRISE  (SO) ,KSTAR (SO) . L. E ,LL,B  (100)  .N.NP.S (SO) . 

1  SURS  (SO)  ,  ALL, INITAL, ABSORB, SIZE,OLD,SBEO,TT,  ASZSI 

REAL***  A  (SOO)  ,  LAN,  V  (40000)  ,  H 

REAL'S  COEPP(IO) ,COV (10,10)  ,U(100)  .PBIRE(IO)  ,IBAR  (10) . 

1  P(100)  .PRVAR(IOO)  , PRBAR  ( 100)  .CORBAR  (100)  .CURVAR  (100) 

C  •••  INITIALIZE  VALUES. 

BEAD  (1,1000)  N, L,K ,LAfl,M,l, ABSORB, 1BITAL, SIZE, TT 
1000  PORBAT  (IS,/,I3./,I6,/,P<*.  1./.P4.  1,/,  I5,/,I5./,IS,/,I5,/,I5) 

NRITE (I, 100S)  N,L,K,LAH,H.X, ABSORB, INITAL, Sill, TT 
1005  PORBAT  (•  N** , IS, •  L«*,I3,*  K**,I6, 

1  •  LAS««,r4.1,*  N*  *  ,  P4.  1,  *  I«*.I5,*  ABSORB**. IS, 

2  •  INITAL**, 15,*  SIZE**, IS,*  TT**,I5) 

NP  *  N  p  1 

BEAD  (11,1020)  SEED 
1020  POBfiAT  (1 10} 

C  •••  FOB  EACH  STATE  DETERBINE  THE  RUBBER  OP  STATES  TBAT  CAN  BE  ENTBBED. 
S(1)  *  1 
DO  BO  J*2,NP 

80  S(J)  =  2 

N8IIE  (3,  10  IS)  (S  ( J),J*1,NP) 

1015  PORBAT  (•  S(J)  **.2015) 

SUBS(I)  *  0 

ALL  *  S(1) 

DO  90  J*2,NP 

SUBS(J)  *  ALL 

80  ALL  *  ALL  P  S(J) 

C  •••  TUE  NEXT  INITIALIZATIONS  BAY  BE  BON  SPECIFIC. 

OLD  *  0 
DO  100  J*2,BP 

C  •••  COBPUTE  TRANSITION  PROBABILITIES. 

P  (SURS  (J)  p  1)  *  N/iLAfltN) 

P  (SUBS  (J)  p2)  *  LAB/  (LAHpN) 

100  CONTINUE 

P(1|  *  1.0 

HRITE  (3,  10  IB)  (P  (LB)  ,  LN*  1,  ALL) 

1016  PORBAT  (•  P  (•)  *  •  ,20FS.  2) 

C  •••  FOR  EACH  STATE  J  DETERBINE  STATES  TBAT  CAB  BE  EBTEBEO. 

DO  400  J*  2,  N 

R  (SURS  (J)  fl)  *  J-2 

4  00  H  (SURS  (J)  p2)  *  J 

B(D  -  1 

B (SUBS  (NP)  p  1)  -  N-1 
fl(SURS  (NP)  p2)  *  N 
NRITE  (3,1017)  (H (LL)  ,LL*  1,  ALL) 

1017  PORBAT  (•  H  (*)  **,20IS) 

L  •  ••  COBPUTE  REUSE  VECTOR. 

00  SOO  J*  1, NP 

A  (ALL*2pSUBS  (J)p  1)  *0.0 


00000100 
00000110 
00000120 
00000130 
00000140 
00000 ISO 
00000160 
00000170 
00000160 
00000190 
00000200 
00000210 
00000220 
00000230 
00000240 
00000250 
00000260 
.00000270 
00000280 
00000290 
00000300 
00000310 
00000320 
00000330 
00000340 
00000 3SO 
00000360 
00000370 
000003B0 
00000390 
00000400 
00000410 
00000420 
00000430 
00000440 
00000450 
00000460 
00000470 
000004B0 
00000490 

ooooosoo 

00000510 
00000S20 
00000530 
00000 S40 
00000550 

ooooosoo 

00000570 

ooooosao 

00000S90 


Fig.  3  Driver  Routine  for  Example 
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500 


1014 

550 

«■  ••• 


A  l ALL* 24 SUBS  (J)  +2) 

A  (AlL*3+SUBS(J)4l| 

A  (ALL*  J+iU  B5  (J|+3) 

a  (sons 

A  (SUBS  (J)42) 

A  (ALL+SUAS 
A  (ALL+SOHS  (JJ  +  2J 
A  (ALLfSOHS  (2)41) 
A(ALL*2+1) 

A  (ALL*24SUBS  (3)  +  1) 
A(ALL*J4SUHS(2)+2) 

A  ( ALL*  J45UBS  (4)  4 1) 

DO  550  LL*1.t 

8B1TE  (3,10181  (A  (ALL* 

rOBAAT  (*  A  (* ,LL)  **, 

CONTINUE 

ASI2E  *  ALL*L 


*0.0 

*0.0 

*0.0 

*  B(SUnS(J)4l)/(tAH48) 

-  B  (SU BS  (J)  42)  /  (LAB 40) 

-  1/(LAB4B) 

*  1/(LAH4V| 

*  1/LAB 

*  1/  (LAB4**) 

*  1/  (LAB+W) 

*  1/(LAB4«) 

-  1/  (LAB40) 


{LL-1)4LBi,LB*1,ALLi 

2QP5.2) 


CALL  PARALLEL  SIBULATION  PBOGSAB. 

CALL  CHAIN  (K.CLD.  1,  NP.IHITAL.  A8S0BB.S,  A.SBBS.ALL.P.ItPBIBE, 

1  KSTAB,li,L.ASI2E,A,PB18E,XBAB,COT.COBPF,TT,PBBAB, 

2  PRVAfi.CUBBAS. CUB VAB, SEED, SIZE, ») 


STOP 

BHD 


00000400 

00000410 

00000420 

00000430 

00000440 

00000450 

00000440 

00000470 

00000400 

00000490 

00000700 

00000710 

00000720 

00000730 

00000740 

00000750 

00000740 

00000770 

00000780 

00000790 

00000800 

00000810 

00000820 


Fig.  3  (Continued) 
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Sample  Output  from  CHAIN  Subroutine 
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«  0.370 *01090-05  0.570595900-05 

2  0.950020000  00  O.0O2*7«35O-OO  0.192775950-00  0. 393030950-00 

3  0.2577712*0  00  0.919551330  00  0.350929000-08  0.707131800-00 
•  0.192905000  00  0.3595*9520  00  0.9710(5010  00  0.151090750-07 
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